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Abstract 

Both orbital and attitude dynamics employ the method of variation of parameters. In 
a non-perturbed setting, the coordinates (or the Euler angles) get expressed as functions 
of the time and six adjustable constants called elements. Under disturbance, each such 
expression becomes ansatz, the “constants” being endowed with time dependence. The 
perturbed velocity (linear or angular) consists of a partial time derivative and a convective 
term containing time derivatives of the “constants.” It can be shown that this construc¬ 
tion leaves one with a freedom to impose three arbitrary conditions upon the “constants” 
and/or their derivatives. Out of convenience, the Lagrange constraint is often imposed. 
It nullifies the convective term and thereby guarantees that under perturbation the func¬ 
tional dependence of the velocity upon the time and “constants” stays the same as in the 
undisturbed case. “Constants” obeying this condition are called osculating elements. 

The “constants” chosen to be canonical, are called Delaunay elements, in the orbital 
case, or Andoyer elements, in the spin case. (As some of the Andoyer elements are time- 
dependent even in the free-spin case, the role of “constants” is played by these elements’ 
initial values.) The Andoyer and Delaunay sets of elements share a feature not readily 
apparent: in certain cases the standard equations render these elements non-osculating. 

In orbital mechanics, elements calculated via the standard planetary equations come 
out non-osculating when perturbations depend on velocities. To keep elements osculating 
under such perturbations, the equations must be amended with extra terms that are not 
parts of the disturbing function (Efroimsky and Goldreich 2003, 2004). For the Kepler el¬ 
ements, this merely complicates the equations. In the case of Delaunay parameterisation, 
these extra terms not only complicate the equations, but also destroy their canonicity. 
So under velocity-dependent disturbances, osculation and canonicity are incompatible. 

Similarly, in spin dynamics the Andoyer elements come out non-osculating under 
angular-velocity-dependent perturbation (a switch to a noninertial frame being one such 
case). Amendment of the dynamical equations only with extra terms in the Hamiltonian 
makes the equations render nonosculating Andoyer elements. To make them osculating, 
more terms must enter the equations (and the equations will no longer be canonical). 

It is often convenient to deliberately deviate from osculation by substituting the La¬ 
grange constraint with an arbitrary condition that gives birth to a family of nonosculat¬ 
ing elements. The freedom in choosing this condition is analogous to the gauge freedom. 
Calculations in nonosculating variables are mathematically valid and sometimes highly 
advantageous, but their physical interpretation is nontrivial. For example, nonosculat¬ 
ing orbital elements parameterise instantaneous conics not tangent to the orbit, so the 
nonosculating inclination will be different from the real inclination of the physical orbit. 

We present examples of situations in which ignoring of the gauge freedom (and of the 
unwanted loss of osculation) leads to oversights. 
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1 Introduction 


1.1 Historical Prelude 

The orbital dynamics is based on the variation-of-parameters method, invention whereof is 
attributed to Euler (1748, 1753) and Lagrange (1778, 1783, 1808, 1809, 1810). Though both 
greatly contributed to this approach, its initial sketch was offered circa 1687 by Newton in his 
unpublished Portsmouth Papers. Very succinctly, Newton brought up this issue also in Cor. 3 
and 4 of Prop. 17 in the first book of his Principia. 

Geometrically, the part and parcel of this method is representation of an orbit as a set of 
points, each of which is contributed by a member of some chosen family of curves C{k), where 
K stands for a set of constants that number a particular curve within the family. (For example, 
a set of three constants k = {a, b, c} defines one particular hyperbola y = ax^ + hx + c 
out of many). This situation is depicted on Fig.l. Point A of the orbit coincides with some 
point Ai on a curve C{ki). Point B of the orbit coincides with point A 2 on some other curve 
C{k 2 ) of the same family, etc. This way, orbital motion from A io B becomes a superposition 
of motion along Ck from Ai to A 2 and a gradual distortion of the curve Ck from the shape 
C{ki) to the shape C{k 2 ). In a loose language, the motion along the orbit consists of steps 
along an instantaneous curve C{k) which itself is evolving while those steps are being made. 



Fig.l. Each point of the orbit is contributed by a member of some family of 
curves C{k) of a certain type, n standing for a set of constants that number 
a particular curve within the family. Motion from A to B is, first, due to the 
motion along the curve C{k) from Ai to A 2 and, second, due to the fact 
that during this motion the curve itself was evolving from C{ki) to C{k 2 ) ■ 
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Normally, the family of curves is chosen to be that of ellipses or that of hyperbolae, k being 
six orbital elements, and A being the time. However, if we disembody this idea of its customary 
implementation, we shall see that it is of a far more general nature and contains three aspects: 

1. A trajectory may be assembled of points contributed by a family of curves of an essen¬ 
tially arbitrary type, not just conics. 

2. It is not necessary to choose the family of curves tangent to the orbit. As we shall see 
below, it is often benehcial to choose them nontangent. We shall also see examples when in 
orbital calculations this loss of tangentiality (loss of osculation) takes place and goes unnoticed. 

3. The approach is general and can be applied, for example, to Euler’s angles. A disturbed 
rotation can be thought of as a series of steps (small turns) along different Eulerian cones. An 
Eulerian cone is an orbit (on the Euler angles’ manifold) corresponding to an unperturbed spin 
state. Just as a transition from one instantaneous Keplerian conic to another is caused by dis¬ 
turbing forces, so a transition from one instantaneous Eulerian cone to another is dictated by 
external torques or other perturbations. Thus, in the attitude mechanics, the Eulerian cones 
play the same role as the Keplerian conics do in the orbital dynamics. Most importantly, a 
perturbed rotation may be “assembled” of the Eulerian cones in an osculating or in a nonoscu¬ 
lating manner. An unwanted loss of osculation in attitude mechanics happens in the same 
way as in the theory of orbits, but is much harder to notice. On the other hand, a deliberate 
choice of nonosculating rotational elements in attitude mechanics may sometimes be benehcial. 

From the viewpoint of calculus, the concept of variation of parameters looks as follows. We 
have a system of differential equations to solve (“system in question”) and a system of differ¬ 
ential equations (“hducial system”) whose solution is known and contains arbitrary constants. 
We then use the known solution to the hducial system as an ansatz for solving the system in 
question. The constants entering this ansatz are endowed with time dependence of their own, 
and the subsequent substitution of this known solution into the system in question yields equa¬ 
tions for the “constants.” The number of “constants” often exceeds that of equations in the 
system to solve. In this case we impose, by hand, arbitrary constraints upon the “constants.” 
For example, in the case of a reduced iV-body problem, we begin with 3(iV — 1) unconstrained 
second-order equations for 3(iV — 1) Cartesian coordinates. After a change of variables from 
the Cartesian coordinates to the orbital parameters, we end up with 3(iV — 1) diherential 
equations for the 6(iV — 1) orbital variables. Evidently, 3(iV — 1) constraints are necessary.^ 
To this end, the so-called Lagrange constraint (the condition of the instantaneous conics being 

^ In a fixed Cartesian frame, any solution to the unperturbed reduced 2-body problem can be written as 
Xj = fj{t, Cl, Ce) , j = 1, 2, 3 , 

Xj = gj{t, Cl, .. .,Cq) , gj = 

the adjustable constants C standing for orbital elements. Under disturbance, the solution is sought as 
^3 = /i(C Ci{t), . . . , Ce{t)) , j = 1, 2, 3 , 

X, = g,it, Ciit), ..., Ceit)) + ^.{t, Ciit), ..., C^it)) , g, = ^ . ^3 = ^ Cr . 

Insertion of xj = fj{t, C) into the perturbed gravity law yields three scalar equations for six functions Cr(t). 
This necessitates imposition of three conditions upon Cr and Cr- Under the simplest choice = 0 , j = 1, 2,3, 
the perturbed physical velocity Xj(t, C) has the same functional form as the unperturbed gj{t, C). Therefore, 
the instantaneous conics become tangent to the orbit (and the orbital elements Cr{t) are called osculating). 
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tangent to the physical orbit) is introduced almost by default, because it is regarded natural. 
Two things should be mentioned in this regard: 

First, what seems natural is not always optimal. The freedom of choice of the supplementary 
condition (the gauge freedom) gives birth to an internal symmetry (the gauge symmetry) of 
the problem. Most importantly, it can be exploited for simplifying the equations of motion for 
the “constants.” On this issue we shall dwell in the current paper. 

Second, the entire scheme may, in principle, be reversed and used to solve systems of differ¬ 
ential equations with constraints. Suppose we have N + M variables Cj{t) obeying a system 
of N differential equations of the second order and M constraints expressed with first-order 
differential equations or with algebraic expressions. One possible approach to solving this sys¬ 
tem will be to assume that the variables Cj come about as constants emerging in a solution to 
some fiducial system of differential equations. Then our N second-order differential equations 
for Cj{t) will be interpreted as a result of substitution of such an ansatz into the hducial 
system with some perturbation, while our M constraints will be interpreted as weeding out of 
the redundant degrees of freedom. Unfortunately, this subject is out of the scope of our paper, 
and it will be developed somewhere else. 


1.2 The simplest example of gauge freedom. 

Variation of constants hrst emerged in the nonlinear context of celestial mechanics and later 
became a universal tool. We begin with a simple example offered in Newman & Efroimsky 
(2003) 

A harmonic oscillator disturbed by a force AF(t) gives birth to the initial-condition problem 
X + X = AF(t) , with a;(0) and i;(0) known , (1) 

whose solution may be sought using ansatz 

X = Ci{t) sint + C 2 {t) cost . (2) 

This will lead us to 

x= Ci{t) sint + C 2 it) cost + Ci{t) cost — C 2 {t) sint . (3) 

It is common, at this point, to put the sum Ci(t) sint -|- C 2 (t) cost equal to zero, in order to 
remove the ambiguity stemming from the fact that we have only one equation for two variables. 
Imposition of this constraint is convenient but not obligatory. A more general way of hxing 
the ambiguity may be expressed as 

C'i(t) sint-|-C' 2 (t) cost = ^(t) , (4) 

0(t) being an arbitrary function of time. This entails: 

X = (j) + C'i(t) cost — (72(t) sint — Cl(t) sint — C 2 (t) cost , (5) 

summation whereof with (2) gives: 

X + X = (j) + Cl(t) cost — C 2 (t) sint . (6) 
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Substitution thereof into (1) yields the dynamical equation re-written in terms of the “con¬ 
stants” Cl, C 2 . This equation, together with identity (4), will constitute the following system: 


This leads to 


0 -|- Ci(t) cost — C 2 (t) sint = AF(t) , 
Ci(t) sint-|-C 2 (t) cost = ^(t) , 

Cl = AF cost — -f- (0cost) 
at 


C 2 = — AF sin t -|- — (0 sin t) , 

(Jj L 

the function 0(t) still remaining arbitrary.^ Integration of (8) entails: 

rt 


Cl = 


AF cost dt 


'Cost -|- ai 


(7) 


( 8 ) 


(9) 


C 2 = — J AF sin t' dt' -1- 0 sin t -|- 02 . 

Substitution of (9) into (2) leads to complete cancellation of the 0 terms: 

a; = Cl sin t -|- C 2 cos t = — cos t J AF sin t' dt' -|- sin t J AF cos t' dt' + Oi sin t -|- 02 cos t (10) 


Naturally, the physical trajectory x{t) remains invariant under the choice of gauge function 
(j){t) , even though the mathematical description (9) of this motion in terms of the parameters 
C is gauge dependent. It is, however, crucial that a numerical solution of the system (8) will 
come out ^-dependent, because the numerical error will be sensitive to the choice of 4>{t). This 
issue is now being studied by P. Gurhl and 1. Klein, and the results are to be published soon. 
(Gurhl & Klein 2006) 

It remains to notice that (8) is a simple analogue to the Lagrange-type system of planetary 
equations, system that, too, admits gauge freedom. (See subsection 2.2 below.) 


1.3 Gauge freedom under a variation of the Lagrangian 

The above example permits an evident extension. (Efroimsky 2002a,b.) Suppose some me¬ 
chanical system obeys the equation 

r = F{t, r, r) , (11) 

whose solution is known and has a functional form 

_ r = fi t, Cl, ..., Ce) , (12) 

^ Function d(t) can afford being arbitrary, no matter what the initial conditions are to be. Indeed, for fixed 
a;(0) and i:(0), the system ^ 2 ( 0 ) = x(0) , ^(0) -I- Gi(0) = x(0) solves for Gi(0) and ^ 2 ( 0 ) for an arbitrary 
choice of d(0). 


Tft€LLhjc.§pL€E.£Lc.€M.L 


5 



Cj being adjustable constants to vary only under disturbance. 

When a perturbation AF gets switched on, the system becomes: 


f = F{t, r, r) + AF{t, r, r) , 
and its solution will be sought in the form of 


Evidently, 


r = fit, C'i(t), . . . , Ceit )) . 


f = A + ^ = yALc- 


( 13 ) 


(14) 


(15) 


In dehance of what the textbooks advise, we do not put # nil. Instead, we proceed further to 


f 


d^f 

dA 


6 

+ E 


i=i 


d^f 

dt dCj 


C, + $ 


(16) 


dot standing for a full time derivative. If we now insert the latter into the perturbed equation 
of motion (13) and if we recall that, according to (11),^ d'^f/dA = F , then we shall obtain 
the equation of motion for the new variables Cj{t) : 


where 


E 


i=i 


d^f A 

dt dCj ' 


+ # 


AF 


$ 


E 




dC, 




(17) 


(18) 


so far is merely an identity. It will become a constraint after we choose a particular functional 

df ■ 

form # (t; Ci, . . . , Cq) for the gauge function # , i.e., if we choose that the sum X) 

be equal to some arbitrarily fixed function #(t; Ci, ..., Cq) of the time and of the variable 
“constants.” This arbitrariness exactly parallels the gauge invariance in electrodynamics: on 
the one hand, the choice of the functional form of #(t; Ci, ..., Cq) will never^ influence the 
eventual solution for the physical variable r ; on the other hand, though, a qualihed choice 
may considerably simplify the process of Ending the solution. To illustrate this, let us denote 
by g(f , Cl, . . . , Ce) the functional dependence of the unperturbed velocity on the time and 
adjustable constants: 


g(i, Cl, .... Ce) = k/(t, Ci, . ... Cs) , (19) 

^ We remind that in (11) there was no difference between a partial and a full time derivative, because at 
that point the integration “constants” Ci were indeed constant. Later, they acquired time dependence, and 
therefore the full time derivative implied in (15 - 16) became different from the partial one implied in (11). 

Our usage of words “arbitrary” and “never” should be limited to the situations where the chosen gauge 
(21) does not contradict the equations of motion (20). This restriction, too, parallels a similar one present in 
field theories. Below we shall encounter a situation where this restriction becomes crucial. 
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and rewrite the above system as 


E 


ds ^ 

dC, 


$ + AF 


( 20 ) 


E 


dC, 


c, 


$ . 


( 21 ) 


If we now dot-multiply the hrst equation with df/dCi and the second one with dg/dCi , 
and then take the difference of the outcomes, we shall arrive at 


j:[CnC,]C, 



df _ ^ ^ 
dC^ ■ dCn 


( 22 ) 


the Lagrange brackets being dehned in a gauge-invariant (i.e., ^-independent) fashion.® If 
we agree that $ is a function of both the time and the parameters Cn, but not of their 
derivatives,® then the right-hand side of (22) will implicitly contain the hrst time derivatives 
of Cn ■ It will then be reasonable to move these to the left-hand side. Hence, (22) will be 
reshaped into 


E 


( 


[Cn C,] 


+ 


dCn 




dCj 

dt 


(9f (9$ (9g 

dCn dCn dt dCn 


(23) 


This is the general form of the gauge-invariant perturbation equations, that follows from 
the variation-of-parameters method applied to problem (13), for an arbitrary perturbation 
F{r, r, t) and under the simplifying assumption that the arbitrary gauge function $ is 
chosen to depend on the time and the parameters Cn , but not on their derivatives.^ As¬ 
sume that our problem (13) is not simply mathematical but is an equation of motion for some 

^ The Lagrange-bracket matrix is defined in a gauge-invariant way: 

rr r ] = 

^ ^ " C - dCn dCj dCn dCj ■ 


and so is its inverse, the matrix composed of the Poisson brackets 


Evidently, (22) yields 


dCn dCj _ dCn ^ 
“ df dg dg df 


Cn = 


■ df 

■ 

(AF - - $ • — 


V J dCj \ 


® The necessity to fix a functional form of $(t; Ci, . . . , Ce) , i.e., to impose three arbitrary conditions 
upon the “constants” Cj , evidently follows from the fact that, on the one hand, in the ansatz (14) we have six 
variables C'„(t) and, on the other hand, the number of scalar equations of motion (i.e., Cartesian projections 
of the perturbed vector equation (13)) is only three. This necessity will become even more mathematically 
transparent after we cast the perturbed equation (13) into the normal form of Cauchy. (See Appendix.) 

^ We may also impart the gauge function with dependence upon the parameters’ time derivatives of all 
orders. This will yield higher-than-first-order derivatives in equation (23). In order to close this system, one 
will then have to impose additional initial conditions, beyond those on r and r . 
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physical setting, so that -F is a physical force corresponding to some undisturbed Lagrangian 
Co 1 and AF is a force perturbation generated by a Lagrangian variation A£. If, for exam¬ 
ple, we begin with Co{r ^ 'f 't t) = ^/2 — U{r, t), momentum p = f, and Hamiltonian 

Ho{r , p, t) = p^/2 -|- U{r , t), then their disturbed counterparts will read: 

r,t) = ^ - U{r) + AC{r,r,t) , (24) 


p = r + 


dAC 


H = p r — C = 


+ u + An , 


An = - AC - 


1 fdAC 


2 V dr 


The Euler-Lagrange equations written for the perturbed Lagrangian (24) are: 




where the disturbing force is given by 


AF = 


dAC d fdAC 


dr dt \ dr 


Its substitution in (23) yields the generic form of the equations in terms of the Lagrangian 
disturbance (Efroimsky & Goldreich 2004): 


E + 


df d fdAC 


dCn dCj \ dr 


d ^ , 1 fdAC 

dCn 2 \ dr 


dg df d d AC d 
dCn dCn dt dr dCn 


This equation not only reveals the convenience of the special gauge 


$ = 


dAC 

dr 


dAC 


(which reduces to $ = 0 in the case of velocity-independent perturbations), but also explicitly 
demonstrates how the Hamiltonian variation comes into play: it is easy to notice that, according 
to (27), the sum in square brackets on the right-hand side of (30) is equal to — An , so the 
above equation takes the form [Cn Cj] Cj = — dAn/dCn ■ All in all, it becomes 

clear that the trivial gauge, # = 0 , leads to the maximal simplihcation of the variation-of- 
parameters equations expressed through the disturbing force: it follows from (22) that 


i:[C„Cj]q = AF 


provided we have chosen # = 0 
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However, the choice of the special gauge (31) entails the maximal simplihcation of the variation- 
of-parameters equations when they are formulated via a variation of the Hamiltonian: 

[Cn Cj] , provided we have chosen $ = — —• (33) 

j diL C/7* 

It remains to spell out the already obvious fact that, in case the unperturbed force F is 
given by the Newton gravity law (i.e., when the undisturbed setting is the reduced two-body 
problem), then the variable “constants” Cn are merely the orbital elements parameterising a 
sequence of instantaneous conics out of which we “assemble” the perturbed trajectory through 
(14). When the conics’ parameterisation is chosen to be via the Kepler or the Delaunay vari¬ 
ables, then (30) yields the gauge-invariant version of the Lagrange-type or the Delaunay-type 
planetary equations, accordingly. Similarly, (22) implements the gauge-invariant generalisation 
of the planetary eqnations in the Euler-Gauss form. 

From (22) we see that the Euler-Gauss-type planetary eqnations will always assnme their 
simplest form (32) nnder the gange choice $ = 0 . In astronomy this choice is called 

“the Lagrange constraint.” It makes the orbital elements osculating, i.e., guarantees that the 
instantaneous conics, parameterised by these elements, are tangent to the perturbed orbit. 

From (33) one can easily notice that the Lagrange- and Delaunay-type planetary eqnations 
simplify maximally nnder the condition (31). This condition coincides with the Lagrange 
constraint $ = 0 when the perturbation depends only upon positions (not npon velocities or 
momenta). Otherwise, condition (31) deviates from that of Lagrange, and the orbital elements 
rendered by equation (33) are no longer osculating (so that the corresponding instantaneous 
conics are no longer tangent to the physical trajectory). 

Of an even greater importance will be the following observation. If we have a velocity- 
dependent pertnrbing force, we can always hnd the appropriate Lagrangian variation and, 
therefrom, the corresponding variation of the Hamiltonian. If now we simply add the negative 
of this Hamiltonian variation to the distnrbing fnnction, then the resulting eqnations (33) will 
render not the osculating elements bnt orbital elements of a different type, ones satisfying 
the non-Lagrange constraint (31). Since the instantaneous conics, parameterised by snch non¬ 
osculating elements, will not be tangent to the orbit, then physical interpretation of snch 
elements may be nontrivial. Besides, they will retnrn a velocity different from the physical 
one.® This pitfall is well camonflaged and is easy to fall in. 

These and other celestial-mechanics applications of the gauge freedom will be considered 
in detail in section 2 below. 


1.4 Canonicity versus osculation 

One more relevant development will come from the theory of canonical perturbations. Snppose 
that in the absence of distnrbances we start ont with a system 


4 


dp 


P 


dq 


(34) 


q and p being the Gartesian or polar coordinates and their conjngated momenta, in the 
orbital case, or the Euler angles and their momenta, in the rotation case. Then we switch, via 

® We mean that substitution of the values of these elements in g(t ; C'i(t), . . . , Ce{t)) will not give the 
right velocity. The correct physical velocity will be given by r = g + $ . 
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a canonical transformation 


q = f{Q, P, t) , p = xiQ, P, t) ( 35 ) 

to 

Q = . P=-_ = 0 . «- = 0. (36) 

where Q and P denote the set of Delaunay elements, in the orbital case, or the initial values 
of the Andoyer variables, in the case of rigid-body rotation. 

This scheme relies on the fact that, for an unperturbed motion (i.e., for an unperturbed 
Keplerian conic, in an orbital case; or for an undisturbed Eulerian cone, in the spin case) a 
six-constant parameterisation may be chosen so that: 

1 . the parameters are constants and, at the same time, are canonical variables {Q, P} 
with a zero Hamiltonian T-C*{Q, P) = 0 ; 

2 . for constant Q and P, the transformation equations (35) are mathematically equiva¬ 
lent to the dynamical equations (34). 


Under perturbation, the “constants” Q, P begin to evolve so that, after their substitution into 
q = f (Qii), Pit), i) , P = xiQii), Pit), i) , (37) 

(/, X being the same functions as in (35)), the resulting motion obeys the disturbed equations 


q = 


d 


P = 


d 


dp dq 

We also want our “constants” Q and P to remain canonical and to obey 


(38) 


Q = 


d {H* + An* 
9P 


P = - 


d {H* + AH*) 


(39) 


where 


P* = 0 and AH* {Q, P t) = An{q{Q,P,t), p{Q,P,t), t) . (40) 

Above all, we wish the perturbed “constants” C = Q, P (the Delaunay elements, in the 
orbital case; or the initial values of the Andoyer elements, in the spin case) to osculate. This 
means that we want the perturbed velocity to be expressed by the same function of Cj(t) and 
t as the unperturbed velocity. Let us check when this is possible. The perturbed velocity is 

q = 9 + ^ (41) 


where 


9{C(t), t) 


dq{C{t), t) 
dt 


(42) 
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is the functional expression for the unperturbed velocity, while 


t) 


E 


t=i 


dq{C{t), t) 

dC, 


Cj{t) 


( 43 ) 


is the convective term. Since we chose the “constants” Cj to make canonical pairs {Q, P) 
obeying (39 - 40), then insertion of (39) into (43) will result in 


^ da ■ ^ 

® = E Qn{t) + E 


9Q 


n=l 


dq 


Pn{t) = 


dAH{q, p) 
dp 


(44) 


So canonicity is incompatible with osculation when AH depends on p. Our desire to keep 
the perturbed equations (39) canonical makes the orbital elements Q , P nonosculating in 
a particular manner prescribed by (44). This breaking of gauge invariance reveals that the 
canonical description is marked with “gauge stiffness” (term suggested by Peter Goldreich). 

We see that, under a momentum-dependent perturbation, we still can use the ansatz (37) for 
calculation of the coordinates and momenta, but can no longer use q = dq/dt for calculating 
the velocities. Instead, we must use q = dq/dt + dAH/dp, and the elements Cj will no 
longer be osculating. In the case of orbital motion (when Cj are the nonosculating Delaunay 
elements), this will mean that the instantaneous ellipses or hyperbolae parameterised by these 
elements will not be tangent to the orbit. (Efroimsky & Goldreich 2003.) In the case of spin, 
the sitnation will be similar, except that, instead of an instantaneous Keplerian conic, one 
will be dealing with an instantaneous Eulerian cone - a set of trajectories on the Euler-angles 
manifold, each of which corresponds to some non-perturbed spin state. (Efroimsky 2004.) 

The main conclusion to be derived from this example is the following: whenever we en- 
connter a distnrbance that depends not only npon positions bnt also upon velocities or mo¬ 
menta, implementation of the afore described canonical-perturbation method necessarily yields 
eqnations that render nonosculating canonical elements. It is possible to keep the elements os¬ 
culating, bnt only at the cost of sacrificing canonicity. For example, nnder velocity-dependent 
orbital pertnrbations (like inertial forces, or atmospheric drag, or relativistic correction) the 
equations for osculating Delannay elements will no longer be Hamiltonian (Efroimsky 2002a,b). 

Above in this snbsection we discussed the disturbed velocity q. How about the disturbed 
momentum? For sufficiently simple nnpertnrbed Hamiltonians, it can be written down very 
easily. For example, for H = Ho + AH = p^/2m + U{q) + AH we get: 


p 


dAC 



9 + ^ + 


dAC 

dq 


9 + 



dAH \ 
dq j 


9 ■ 


(45) 


In this case, the pertnrbed momentnm p coincides with the nnpertnrbed one, g. In application 
to the orbital motion, this means that contact elements (i.e., the nonosculating orbital elements 
obeying (31)), when snbstitnted in g(t ] Ci, , Ge), furnish not the correct pertnrbed velocity 

but the correct perturbed momentum, i.e., they osculate the orbit in phase space. Existence 
of snch elements was pointed out long ago by Goldreich (1965) and Brumberg et al (1971). 
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2 


Gauge freedom in the theory of orbits. 

2.1 Geometrical meaning of the arbitrary gauge function $ 

As explained above, the content of subsection 1.3 becomes merely a formulation of the Lagrange 
theory of orbits, provided F stands for the Newton gravity force, so that the undisturbed 
setting is the two-body problem. Then (22) expresses the gauge-invariant (i.e., taken with 
an arbitrary gauge ; Ci, ..., G)) planetary equations in the Euler-Gauss form. These 
equations render orbital elements that are, generally, not osculating. Equation (32) stands for 
the customary Euler-Gauss-type system for osculating (i.e., obeying # = 0 ) orbital elements. 

Similarly, equation (30) stands for the gauge-invariant Lagrange-type or Delaunay-type 
(dependent upon whether G* stand for the Kepler or Delaunay variables) equations. Such 
equations yield elements, which, generally, are not osculating. In those equations, one could 
hx the gauge by putting $ = 0 , thus making the resulting orbital elements osculating. How¬ 
ever, this would be advantageous only in the case of velocity-independent AC . Otherwise, a 
maximal simplihcation is achieved through a deliberate refusal from osculation: by choosing 
the gauge as (31) one ends up with simple equations (33). Thus, gauge (31) simplihes the 
planetary equations. (See equations (46 - 57) below.) Besides, in the case when the Delau¬ 
nay parameterisation is employed, this gauge makes the equations for the Delaunay variables 
canonical for reasons explained above in subsection 1.4. 

The geometrical meaning of the convective term $ becomes evident if we recall that a 
perturbed orbit is assembled of points, each of which is donated by one representative of a 
sequence of conics, as on Fig.2 and Fig.3 where the “walk” over the instantaneous conics may 
be undertaken either in a non-osculating manner or in the osculating manner. The physical 
velocity r is always tangent to the perturbed orbit, while the unperturbed Keplerian velocity 
g = df/dt is tangent to the instantaneous conic. Their difference is the convective term #. 
So if we use non-osculating orbital elements, then insertion of their values in f(t ; Gi, ..., C^) 
will yield a correct position of the body. However, their insertion in g(t; Gi, ..., Cq) will 
NOT give the right velocity. To get the correct velocity, one will have to add $ . (See Appendix 
for a more formal mathematical treatment in the normal form of Gauchy.) 

When using non-osculating orbital elements, we must always be careful about their physical 
interpretation. On Fig. 2, the instantaneous conics are not supposed to be tangent to the 
orbit, nor are they supposed to be even coplanar thereto. (They may be even perpendicular 
to the orbit! - why not?) This means that, for example, the non-osculating element i may 
considerably differ from the real, physical inclination of the orbit. 

We would add that the arbitrariness of choice of the function ^>(t, Gi(t), ..., C^it)) had 
been long known but never used in astronomy until a recent effort undertaken by several authors 
(Efroimsky 2002a,b ; Newman & Efroimsky 2003; Slabinski 2003; Efroimsky & Goldreich 
2003, 2004; Gurhl 2004; Efroimsky 2005c) Substitution of the Lagrange constraint $ = 0 
with alternative choices does not influence the physical motion, but alters its mathematical 
description (i.e., renders different values of the orbital parameters Ci{t)). Such invariance of a 
physical picture under a change of parameterisation goes under the name of gauge freedom. It 
is a part and parcel of electrodynamics and other held theories. In mathematics, it is described 
in terms of hber bundles. A clever choice of gauge often simplihes solution of the equations 
of motion. On the other hand, the gauge invariance may have implications upon numerical 
procedures. We mean the so-called ’’gauge drift,” i.e., unwanted displacement in the gauge 
function $, caused by accumulation of numerical errors in the constants. 
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Fig.2. The orbit is a set of points, each of which is donated by one of the 
confocal instantaneous ellipses that are not supposed to be tangent or even 
coplanar to the orbit. As a result, the physical velocity r (tangent to the orbit) 
differs from the unperturbed Keplerian velocity g (tangent to the ellipse). To 
parameterise the depicted sequence of non-osculating ellipses, and to single it 
out of the other sequences, it is suitable to employ the difference between r and 
g, expressed as a function of time and six (non-osculating) orbital elements: 
$(t, Cl, ..., Ce) = r(t, Cl, ..., Ce) - g(t , Cl, ..., Ce) . Evidently, 


r 


^ _L "V C 
j=i 


g + $ 


where the unperturbed Keplerian velocity is g = drjdt. The convective term, 
which emerges under perturbation, is $ = X) [dr/dCj) Cj. When a particular 
functional dependence of $ on time and the elements is fixed, this function, 
$(t , Cl , ..., Cg), is called gauge function or gauge velocity or, simply, gauge. 



Fig.3. The orbit is represented by a sequence of confocal instantaneous ellipses 
that are tangent to the orbit, i.e., osculating. Now, the physical velocity r 
(tangent to the orbit) coincides with the unperturbed Keplerian velocity g 
(tangent to the ellipse), so that their difference $ vanishes everywhere: 

$(t, Ci,...,C6)= r(t,Ci,...,C6)-g(t,Ci,...,C6)= 5]^C, =0. 

This equality, called Lagrange constraint or Lagrange gauge, is the necessary 
and sufficient condition of osculation. 
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2.2 Gauge-invariant planetary equations 
of the Lagrange and Delaunay types 

We present the gauge-invariant Lagrange- and Delaunay-type equations, following Efroimsky & 
Goldreich (2003). These equations follow from (30) if we take into account the gauge-invariance 
(i.e., the ^-independence) of the Lagrange-bracket matrix [CiCj]. 


da 

dt 


2_ p(-A7f) d^C d / d^C \ 
na dMo dr dMo \ (9f J 



aA£ \ 
dr J 


dg 

dMo 


hh d A + 

dMo dt \ dr J 


(46) 


de 1 — d{ — AH) 
dt na^ e dMo 


dAC d / dAC \ 
dr (9a y dr J 



dAC 

dr 


(1 - e^)V^ 
na‘^ e 


A_ d A + 

dMo dMo dt \ 

d (- AH) _ dAC d 
duj dr duj 


dAC \ 
dr j 



dAC \ 

dr J 



dAC \ 

dr J 


doj 


A d A + ^)1 

dt\ df l\ ' 


(47) 


duj 

dt 


— cos i 

na^ (1 — sin i 


d (- AH) 
di 



dAC 

dr 


dAC ^ / dAC \ 
dr (9i y ^ dr j 


di 


A d A + 

di dt\ dr J 


+ 


(1 - e^)V^ \ d{-AH) 
e de 


dAC ^ / ,9A£ \ 

dr (9e y ^ dr ) 


(48) 



,9A£ \ 

dr J 


de 


A d A + ddd) 

de dt \ dr J 
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di cos i 

dt na^{l — sin i 


dAC\ dg 


na^ (1 


dr / doj 


I _ I d {-AH) 

e2)i/2sini d^l 


na^ (1 — sin i 


dr / di 


dAC 

d 

dr 

duj 

^ . 

d / 

doj ( 

it y 

dAC 

d 

dr 

dn 

dg 

Of 

dn 

dVi 

dAC 

d 

dr 

di 

_ ^ 

d 

di 

dt 


dAC 
^ dr 

dAC\ 


dAC 

dr 


dAC 

dr 


dAC 

dr 


dAC 

dr 


1 - \d{-Al-L) dAC d 


no? e 


dr de 


dAC 


dr / de de dt I dr 


yi_ I d {-AH) _ dAC d 
na da dr da 


dAC 


* + 1 fs + 

dr / da da dt 1 dr 


T^uJmjc. Tft€LLhjc.§pL€E.£LC.€M.L ’^ti.y.A.tc.A. 



Similarly, the gauge-invariant Delaunay-type system can be written down as: 


dL d{ — AH) dAL d / (9A£\ / dAC\ dg dr d / dAC\ 

dt dMo dr dMo \ (9r j (9r j dMo dMo dt \ dr ) 

dMo d{ — AH) dAC d / dAC\ dAC\ dg dr d / dAC\ 

dt dL ^ dr dL \ (9r J \ ^ dr j dL ^ dL dt \ dr ) 

^ _ d (- AH) _ dAC / dAC \ _ / dAC \ ^ ^ ^ dAC \ 

dt doj dr dijj\ ^ dr j y (9f j doj doj dt\ dr j 

doj d{ — AH) dAC d / dAC\ / dAC\ dg dr d / dAC\ 

H ^ ^ ^ ^ j ^ ^ dr j'^ ^ '^Jt[ ^ dr ) 

m _ d (- AH) _ dAC / dAC \ _ / dAC \ ^ dAC \ 

dt dVL dr dVL y (9r J y (9r J <9^ dVL dt y dr ) 

dVt d{ — AH) dAC d / dAC\ / dAC\ dg dr d / dAC\ 
dt dH dr dH y ^(9ry~*~y ^ dr j dH dH dt y ctr J 

where p stands for the reduced mass, while 

L = y-'=o‘/", G = (l - , H = {l - e^y'" cost . (58) 


(52) 

(53) 

. (54) 
, (55) 

(56) 
^ (57) 


The symbols /, g now denote the functional dependencies of the gauge, position, and 
velocity upon the Delaunay, not Keplerian elements, and therefore these are functions different 
from #, /, g used in (46 - 51) where they stood for the dependencies upon the Kepler ele¬ 
ments. (In Efroimsky (2002a,b) the dependencies #, /, g upon the Delaunay variables were 
equipped with tilde, to distinguish them from the dependencies upon the Kepler coordinates.) 

To employ the gauge-invariant eqnations in analytical calculations is a delicate task: one 
should always keep in mind that, in case # is chosen to depend not only npon time bnt also 
npon the “constants” (but not upon their derivatives), the right-hand sides of these eqnations 
will implicitly contain the first derivatives dCi/dt , and one will have to move them to the left- 
hand sides (like in the transition from (22) to (23)). The choices $ = 0 and# = —dAC/dr 
are exceptions. (The most general exceptional gauge reads as # = — dAC/dr + rjit) , 
where r]{t) is an arbitrary function of time.) 

As was expected from (30), both the Lagrange and Delannay systems simplify in the gauge 
(31). Since for orbital motions we have dLL/dp = — dAC/dr , then (31) coincides with (44). 
Hence, the Hamiltonian analysis (34 - 44) explains why it is exactly in the gange (31) that the 
Delannay system becomes symplectic. In physicists’ parlance, the canonicity condition breaks 
the gauge symmetry by stiffly fixing the gange (44), gange that is eqnivalent, in the orbital 
case, to (31) - phenomenon called “gange stiffness.” The phenomenon may be looked npon 
also from a different angle. Above we emphasised that the gange freedom implies essential 
arbitrariness in our choice of the functional form of #(t; Ci, . . . , C^) , provided the choice 
does not come into a contradiction with the equations of motion - an important clause that 
shows its relevance in (34 - 44) and (51 - 56): we see that, for example, the Lagrange choice 
# = 0 (as well as any other choice different from (31)) is incompatible with the canonical 
structure of the equations of motion for the elements. 
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3 A practical example on gauges: 

a satellite orbiting a processing oblate planet. 

Above we presented the Lagrange- and Delaunay-type planetary equations in the gauge- 
invariant form (i.e., for an arbitrary choice of the gauge function ; Ci, Cq)) and 
for a generic perturbation AC that may depend not only upon positions but also upon veloc¬ 
ities and the time. We saw that the disturbing function is the negative Hamiltonian variation 
(which differs from the Lagrangian variation when the perturbation depends on velocities). 
Below, we shall also see that the functional dependence of ATi. upon the orbital elements is 
gauge-dependent. 

3.1 The Gauge Freedom and the Freedom of Frame Choice 

In the most compressed form, implementation of the variation-of-constants method in orbital 
mechanics looks like this. A generic solution to the two-body-problem is expressed with 

r = (59) 


[dtJc P f 

and is used as an ansatz to describe the perturbed motion: 
r = t) , 


dCi dt 


dg 


dg dCj 
dCi dt 


P f 


dg dCj 

dCi dt 


dt 


As can be seen from (63), our choice of a particular gauge is equivalent to a particular way of 
decomposition of the physical motion into a movement with velocity g along the instantaneous 
conic, and a movement caused by the conic’s deformation at the rate $. Beside the fact that 
we decouple the physical velocity r in a certain proportion between these two movements, 
g and $, it also matters what physical velocity (i.e., velocity relative to what frame) is 
decoupled in this proportion. Thus, the choice of gauge does not exhaust all freedom: one can 
still choose in what frame to write ansatz (62) - one can write it in inertial axes or in some 
accelerated or/and rotating ones. For example, in the case of a satellite orbiting a processing 
oblate primary it is most convenient to write the ansatz for the planet-related position vector. 

The kinematic formulae (62 - 64) do not yet contain information about our choice of the 
reference system wherein to employ variation of constants. This information shows up only 
when (62) and (64) get inserted into the equation of motion r + (/rr/r^) = AF to render 


dg dCj 
dCi dt 


= Af = 


dAC 

dr 


dAC 

dr 
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Information about the reference frame, where we employ the method and dehne the elements 
Ci , is contained in the expression for the perturbing force A/. If the operation is carried 
out in an inertial system, A/ contains only physical forces. If we work in a frame moving 
with a linear acceleration a, then A/ also contains the inertial force —a. In case this 
coordinate frame also rotates relative to inertial ones at a rate /x, then A/ also includes the 
inertial contributions —2/x x r — fi x r — x {/j, x r). When studying orbits about an 
oblate precessing planet, it is most convenient (though not obligatory) to apply the variation- 
of-parameters method in axes coprecessing with the planet’s equator of date: it is in this 
coordinate system that one should write ansatz (62) and decompose r into g and $. This 
convenient choice of coordinate system will still leave one with the freedom of gauge nomination; 
in the said coordinate system, one will still have to decide what function $ to insert in (63). 


3.2 The disturbing function in a frame 

co-precessing with the equator of date 

The equation of motion in the inertial frame is 

r" = 

dr ’ 

where U is the total gravitational potential, and time derivatives in the inertial axes are denoted 
by primes. In a coordinate system precessing at angular rate ^{t) , equation (66) becomes: 



.. dU . . , , 

r = —— — 2u, X r — u X r — u, X [u X r) 

dr 


dr 


dAU 

dr 


— 2iJ, X r — fi X r — X X r) 


(67) 


dots standing for time derivatives in the co-precessing frame, and n being the coordinate 
system’s angular velocity relative to an inertial frame. Formula (125) in the Appendix gives 
the expression for /x in terms of the longitude of the node and the inclination of the equator of 
date relative to that of epoch. The physical (i.e., not associated with inertial forces) potential 
U{r) consists of the (reduced) two-body part Uoir) = —GM r/r^ and a term AU{r) 
caused by the planet’s oblateness (or, generally, by its triaxiality). 

To implement variation of the orbital elements dehned in this frame, we note that the 


disturbing force on 

the right-hand side of (67) is generated, according to (65), by 


AC 

(r, r, t) = — AU{r) + r {^ x r) + ^ (^t x r)-(^t x r) . 

(68) 

Since 

dAC 

(69) 


n- = ^ X r , 

ar 

then 

dAC 

(70) 


p = r + = r + p X r 

dr 
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and, therefore, the corresponding Hamiltonian perturbation reads: 


1 /dAC 

2 I dr 


AU + p ■ X r) ] = 


AU + {r X p) ■ p,] = AU — J ■ p , 


with vector J = r x p being the satellite’s orbital angular momentum in the inertial frame. 
According to (63) and (70), the momentum can be written as 

p = g + ^ + pxf , (72) 

whence the Hamiltonian perturbation becomes 


AH = - AC 


1 fdAC 


2 \ dr 


AU + {f xg) ■ p + {^ + px f) ■ {px f)] . (73) 


This is what one is supposed to plug in (30) or, the same, in (46 - 57). 


3.3 Planetary equations in a precessing frame, 
written in terms of contact elements 

In the preceding subsection we fixed our choice of the frame wherein to describe the orbit. 
By writing the Lagrangian and Hamiltonian variations as (68) and (73), we stated that our 
elements would be dehned in the frame coprecessing with the equator. The frame being fixed, 
we are still left with the freedom of gauge choice. As evident from (33) or (46 - 57), the special 
gauge (31) ideally simplihes the planetary equations. Indeed, (31) and (69) together yield 


$ = 


dAC 


= -pxr = -px f , 


wherefrom the Hamiltonian (73) becomes 


= - [ - AU{f) + /X ■ (/ X g) 


while the planetary equations (30) get the shape 


[Cr Q] 


da _ <9 ( - ) 


or, the same. 


[C, Q] — = -^ 
^ dt dCr 


AU{f ) + p -if xg) 


where f and g stand for the undisturbed (two-body) functional expressions (59) and (60) 
of the position and velocity via the time and the chosen set of orbital elements. Planetary 
equations (76) were obtained with aid of (74), and therefore they render non-osculating orbital 
elements that are called contact elements. This is why we equipped the Hamiltonian (75) with 
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superscript '^(cont)." In distinction from the osculating elements, the contact ones osculate 
in phase space: (72) and (74) entail that p = g . As already mentioned in the end of 
section 1, existence of such elements was pointed out by Goldreich (1965) and Brumberg et 
al (1971) long before the concept of gauge freedom was introduced in celestial mechanics. 
Brumberg et al (1971) simply defined these elements by the condition that their insertion 
in g(t ] Cl, , Cq) returns not the perturbed velocity, but the perturbed momentum. 

Goldreich (1965) dehned these elements (without calling them “contact”) differently. Having 
in mind inertial forces (67), he wrote down the corresponding Hamiltonian (71) and added 
its negative to the disturbing function of the standard planetary equations (without enriching 
the equations with any other terms). Then he noticed that those equations furnished non¬ 
osculating elements. Now we can easily see that both Goldreich’s and Brumberg’s dehnitions 
follow from the gauge choice (31). 

When one chooses the Keplerian parameterisation, then (77) becomes: 


d ( - 

m'o 


de 1-e^ (9 ( - (i _ e2)V2 9 ( - 

dt na‘^ e OMo na^ e du 


du -cos^ (i_e2)i/2 

dt na‘^(l — e^)G2 sin i di na'^e de 


dt na^ (1 — e2)H2 sin i 


d ( - A7f(" 


na^ (1 — e^)^A sin i 


a (- AH<' 


. ( 81 ) 


(1 — e^)^A sin i 


d ( - A7f(^ 


d ( - A7f(^ 


no? e 


2 ( 9 (- 


The above equations implement an interesting pitfall. When describing orbital motion relative 
to a frame coprecessing with the equator of date, it is tempting to derive the Hamiltonian 
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variation caused by the inertial forces, and to simply plug it, with a negative sign, into the dis¬ 
turbing function. This would entail equations (76 - 83) which, as demonstrated above, belong 
to the non-Lagrange gauge (31). The elements furnished by these equations are nonosculating, 
so that the conics parameterised by these elements are not tangent to the perturbed trajectory. 
For example, i gives the inclination of the instantaneous non-tangent conic, but differs from 
the real, physical physical (i.e., osculating), inclination of the orbit. This approach - when an 
inertial term is simply added to the disturbing function - was employed by Goldreich (1965), 
Brumberg et al (1971), and Kinoshita (1993), and many others. Goldreich and Brumberg 
noticed that this destroyed osculation, Kinoshita failed to. 

Goldreich (1965) studied how the equinoctial precession of Mars influences the long-term 
evolution of Phobos’ and Deimos’ orbit inclinations. Goldreich assumed that the elements 
a and e stay constant; he also substituted the Hamiltonian variation (75) with its orbital 
average, which made his planetary equations render the secular parts of the elements. He 
assumed that the averaged physical term (AG) is only due to the primary’s oblateness: 


(AG) 


J 2 2 3 cos^ i — 1 

~ (1 - 


(84) 


p being the mean equatorial radius of the planet,® and n being the satellite’s mean motion. 
To simplify the inertial term, Goldreich employed the well known formula 


where 


r X g 


= \jGm a (1 — e^) 


w 


(85) 


in = xi sin i sin G — X 2 sin i cos G -|- X 3 cos i 


( 86 ) 


is a unit vector normal to the instantaneous ellipse, expressed through unit vectors xi, X 2 , X 3 
associated with the co-precessing frame xi, X 2 , X 3 (axes xi and X 2 lying in the planet’s 
equatorial plane of date, and xi pointing along the hducial line wherefrom the longitude of 
the ascending node of the satellite orbit, G , is measured). This resulted in 


^^U^cont)) ^ _ [_(AG) + (M-(/Xg))] = 


Gm J 2 pI 3 cos^ i — 1 
4 (^1 _ 



(/ii sin i sin G 


p ,2 sin i cos VL + p,'^ cos i) 


(87) 


all letters now standing not for the appropriate variables but for their orbital averages. Sub¬ 
stitution of this averaged Hamiltonian in (81 - 82) lead Goldreich, in assumption that both 
|/i|/ (J 2 sin i) and |/x|/ (n J 2 sin i) are much less than unity, to the following system: 


dVL 3 ^ fPe\‘^ cos i 

dt 2 ^ V a y (1 - e^f 


( 88 ) 


di 

dt 


K, — Pi cos Vt — P 2 sin G 


(89) 


Goldreich used the nonsphericity parameter J = (3/2) [pe/p)^ J 2 , where pe is the mean equatorial radius. 
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whose solution, 


i = - — cos [ - X (t - to) + nj + — sin [ - X (t - to) + ^o] + io 
X X 


n = - X (t - to) + fio 


where x 



cos i 

(1 - 


(90) 


tells us that in the course of equinoctial precession the satellite inclination oscillates abont io ■ 
Goldreich (1965) noticed that his i and the other elements were not osculating, bnt he 
assnmed that their secular parts would differ from those of the osculating ones only in the 
orders higher than 0(|/x|). Below we shall probe the applicability limits for this assumption. 
(See the end of snbsection 3.5.) 


3.4 Planetary equations in a precessing frame, 
in terms of osculating elements 

When one introdnces elements in the precessing frame and also demands that they oscnlate in 
this frame (i.e., obey the Lagrange constraint = 0), the Hamiltonian variation readsd*^ 

= - I - AC/ + M ■ (/ X g) + (m X /) ■ (m X /) ] , (91) 


while equation (30) becomes: 


[Cn Q] 


dCj 

dt 


M- 


dC„. 


X g - / X 


dg 

dCn 


A- / X 


dAH^osc) 

dCo. 




(92) 


To ease the comparison of this eqnation with (77), it is convenient to split the expression (91) 
for into two parts: 


= - I R^aUf, «) + M ■ (/ X g) ] 


(93) 


and 


- (/l X /) ■ (^ X /) , 


(94) 


and then to gronp the latter part with the last term on the right-hand side of (35): 


[Cn Q] 


dCj 

dt 


dCn 


+ M ■ 



X g - / X 



A- / X 


dCn 


+ (m X /) 


d 


dCr 


(m X /) 


(95) 


10 Both ATILo"*) and are equal to - [ - At/(/, t) + p • J] = - [ - AU{f, t) + p • (/ x p)] ■ 

However, the canonical momentum now is different from g and reads as: p = g + (p x /) . Hence, the 
functional forms of A7i:(°^=)(/, p) and A 7 i:(“») (/^ p) are different, though their values coincide. 
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Comparison of this analytical theory with a straightforward numerical integration^^ has con- 
hrmed that the 0(|/xp) term in (95) may be neglected over time scales of, at least, hundreds 
of millions of years. In this approximation there is no the difference between and 

, so we shall write down the equations as: 


[Cn Q] 


dCj 

dt 


dC„, 


+ M ■ 


dC^. 


X g - / X 


dC„ 


- A- / X 


dCr,. 


(96) 


For Ci chosen as the Kepler elements, inversion of the Lagrange brackets in (90) will yield 
the following Lagrange-type system: 


da 2 

dt na 




cont) 


dMo 


A- / X 


df 

dM„ 


(97) 


de 1 — 


dt na^ e 


d ( - A7f( 


cont) 


DM, 


- A- / X 


df 

dM„ 


(98) 


(1 - 
na^ e 


d ( - A7f(' 


cont) 


doj 


. (df n dg'' 

+ A^- ^xg-/x — 
\OU OUJ . 


• (f 


duj 

dt 


— cos I 


na^{l — sin i 




di 


(df . dg\ . / df 


(99) 


(l-e^)V^ 
na^ e 


d ( - AH^ 


cont) 


de 


, (df » dg\ . / df" 

+ M- ^xg-/x— - n - [f X — 

\ de de \ de , 


di cos i 

dt na^ (1 — sin i 


a(_AH(“”‘>) (8f ^ ds\ . (, df' 

- -b- - + ^‘ o^^S-Zx— -M' /x — 

du \du doj \ duj. 


( 100 ) 


2 (1 — e^y/"^ sin i 


na 


d ( - AH^ 


cont) 


dn 


(df 


^ dn 




Credit for this comparison goes to Pini Gurfil and Valery Lainey. 
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dn _ 1 [ <9 (- 

dt na^ (1 — sin i di 




1 - 

na^ e 


d ( - 


d ( - A7^(" 


(df dg\ . / df 

+ ^Xg-/X- -M' /X^ 


( 102 ) 


(df dg\ .( df 

\ oa oa \ oa 


terms ^t■( {df/dMo) x g - {dg/dMo) X / ) being omitted in (97 - 98), because these terms 
vanish identically (see the Appendix to Efroimsky (2004)). 


3.5 Comparison of calculations performed in the two above gauges 

Simply from looking at (76 - 83) and (96 - 102) we notice that the difference in orbit descriptions 
performed in the two gauges emerges already in the first order of the precession rate and 
in the hrst order of /7 . 

Calculation of the fi- and /x-dependent terms emerging in (97 - 102) takes more than 20 
pages of algebra. The resulting expressions are published in Efroimsky (2005a), their detailed 
derivation being available in web-archive preprint Efroimsky (2004). As an illustration, we 
present a couple of expressions: 


(1 - ^ 

(1 -1- e cosz/)" 


{ /ii [ — cosh! sin(ci; + u) — sin hi cos(ci; -1- z/) cos i ] sin(ci; -|- u) 


+ f 2 [ — sin n sin(ci; + u) + cos hi cos(a; -|- z/) cos i ] sin(a; -|- z/) 


-|- /is sin(a; -|- u) cos(a; -|- z/) sin i } , 


(103) 


M- ^xg-/x — 
V oe oe 


na^ (3 e -7 2 cos z/ 4- cos u) 
(1 4- e cosz/) Vl — 


(104) 
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V denoting the true anomaly. The fact that almost none of these terms vanish reveals that 
equations (76 - 83) and (96 - 102) may yield very different results, i.e., that the contact elements 
may differ from their osculating counterparts already in the hrst order of /x . 

Luckily, in the practical situations we need not the elements per se but their secular parts. 
To calculate these, one can substitute both the Hamiltonian variation and the /x- and /x- 
dependent terms with their orbital averages^^ calculated through 

(1 - /•2- dv 

2 TT Xo (1 + e cos v'f 



The situation might simplify very considerably if we could also assume that the precession rate 
/X stays constant. Then in equations (97 - 102), we would assume /x = const and proceed 
with averaging the expressions ( {dfjdCj) x g — / x [dg/dCj) ) only (while all the terms 
with /X will now vanish). 

Averaging of the said terms is lengthy and is presented in the Appendix to Efroimsky 
(2004). All in all, we get, for constant /x ; 





g - / 



df ^ . dg\ 3 

& as =2 




IG m (1 — e2) 


,(106) 


^ ^ ^ X g - / X ^ ^ ) = 0 , Cj = e, n, u , i, Mo . (107) 

Since the orbital averages (107) vanish, then e will, along with a , stay constant for as long 
as our approximation remains valid. Besides, no trace of /x will be left in the equations for 
H and i . This means that, in the assumed approximation and under the extra assumption 
of constant /x , the afore quoted analysis (84 - 90), offered by Goldreich (1965), will remain 
valid at time scales which are not too long. 

In the realistic case of time-dependent precession, the averages of terms containing /x 
and /X do not vanish (except for /x- ( {df /dMo) x g — / x (dg/OMo) ) , which is identically 
nil). These terms show up in all equations (except in that for a ) and influence the motion. 
Integration that includes these terms gives results very close to the Goldreich approximation 
(approximation (90) that neglects the said terms and approximates the secular parts of the 
nonosculating elements with those of their osculating counterparts). However, this agreement 
takes place only at time scales of order millions to dozens of millions of years. At larger time 
scales, differences begin to accumulate (Lainey et al 2005). 

In real life, the equinoctial-precession rate of the planet, /x , is not constant. Since the 
equinoctial precession is caused by the solar torque acting on the oblate planet, this precession 
is regulated by the relative location and orientation of the Sun and the planetary equator. This 
is why /X of a planet depends upon this planet’s orbit precession caused by the pull from the 
other planets. This dependence is described by a simple model developed by Golombo (1966). 

^2 Mathematically, this procedure is, to say the least, not rigorous. In practical calculations it works well, at 
least over not too long time scales. 
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4 Conclusions: how we benefit from the gauge freedom. 

In this article we gave a review of the gauge concept in orbital and attitude dynamics. Essen¬ 
tially, this is the freedom of choosing nonosculating orbital (or rotational) elements, i.e., the 
freedom of making them deviate from osculation in a known, prescribed, manner. 

The advantage of elements introduced in a nontrivial gauge is that in certain 
situations the choice of such elements considerably simplifies the mathematical 
description of orbital and attitude problems. One example of such simplihcation is the 
Goldreich (1965) approximation (90) for satellite orbiting a processing oblate planet. Although 
performed in terms of non-osculating elements, Goldreich’s calculation has the advantage of 
mathematical simplicity. Most importantly, later studies (Efroimsky 2005a,b) have conhrmed 
that Goldreich’s results, obtained for nonosculating elements, serves as a very good approxi¬ 
mation for the osculating elements. To be more exact, the secular parts of these nonosculating 
elements coincide, in the hrst order over the precession-caused perturbation, with those of their 
osculating counterparts, the difference accumulating only at very long time scales - see the end 
of section 3 above. A comprehensive investigation into this topic, with the relevant numerics, 
will be presented in Lainey et al (2005). 

On the other hand, neglect of the gauge freedom may sometimes produce camou¬ 
flaged pitfalls caused by the fact that nonosculating elements lack evident physical 
meaning. For example, the nonosculating “inclination” does not coincide with the real, physi¬ 
cal inclination of the orbit. This happens because nonosculating elements parameterise instan¬ 
taneous conics nontangent to the orbit. Similarly, nonosculating Andoyer elements L , G , H 
are no longer the same projections of the angular momentum as their osculating counterparts. 


Appendix 1. Mathematical formalities: 

Orbital dynamics in the normal form of Cauchy 


Let us cast the perturbed equation 



f = 

F + Af - 

= — 

^7 + 

(108) 

into the normal 

form of Gauchy: 





r = 

V , 




(109) 

V = 


( r{t, Cl , ... 


) , v(t, Gi, ..., Ge) ) . 

(110) 

Insertion of our 

ansatz 






r = 

fit, Ci{t), 


, c,it)) , 

(111) 

will make (109) 

equivalent to 







df ^ 

E 

dC, • 

(112) 
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The function / is, by definition, the generic solution to the unperturbed equation 


This circumstance, along with (112), will transform (109) into 

i: ^ C. + <i> = AF( /(«, Ci,...,Ce) , g(«, Ci,...,Ce) + # ) 

i ^ 

where 

(115) 

is an identity, fit, Ci , ..., C^) and g(t, Ci , ..., Cf) = df/dt being known functions. Now 
(114 - 115) make an incomplete system of six hrst-order equations for nine variables {Ci , ... , Cg , 
$ 1 , ... , $ 3 ). So one has to impose three arbitrary conditions on C*, , for example as 

$ = #(f, C'i,...,C'6 ) . (116) 


(113) 

(114) 


This will result in a closed system of six equations for six variables Cj : 


C. = AF( /(«. Ci, ....Cb) , g(«, Ci, ....Cb) + <6 ) - 6 


^ ^ dCj 
, dCi dt 


(117) 

(118) 


$ = ^>( t, Cl, ..., Cg ) now being some hxed function (gauge).A trivial choice is 

$( f, Cl, ..., Cg ) = 0, and this is what is normally taken by default. This choice is only 
one out of inhnitely many, and often is not optimal. Under an arbitrary, nonzero, choice of the 
function $( t, Ci, ..., Cg ), the system (117 - 118) will have some different solution Cj{t) . 
To get the appropriate solution for the Cartesian components of the position and velocity, one 
will have to use formulae 


r = /(t, Ci,...,Cg) 


(119) 


r 


V = g( t, Cl, ..., Cg ) + $( t, Cl, ..., Cg ) , 


( 120 ) 


Generally, $ may depend also upon the variables’ time derivatives of all orders: $(t; Ci, Ci, Ci, ...). 
This will give birth to higher time derivatives of C in subsequent developments and will require additional 
initial conditions, beyond those on r and r , to be fixed to close the system. So it is practical to accept (116). 
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Appendix 2. Precession of the equator of date 


relative to the equator of epoch 

The afore introduced vector /x is the precession rate of the equator of date relative to 
the equator of epoch. Let the inertial axes (X , Y , Z) and the corresponding unit vectors 
(X , Y , Z ) be hxed in space so that X and Y belong to the equator of epoch. A rotation 
within the equator-of-epoch plane by longitude hp , from axis X , will dehne the line of nodes, 
X . A rotation about this line by an inclination angle Ip will give us the planetary equator 
of date. The line of nodes x , along with axis y naturally chosen within the equator-of-date 
plane, and with axis orthogonal to this plane, will constitute the precessing coordinate 
system, with the appropriate basis denoted by (x, y , z ) . 

In the inertial basis (X, Y , Z ) , the direction to the North Pole of date is given by 

T 

z = (sinJpSinhp , — sin/pCoshp , cos Jp) (121) 

while the total angular velocity reads: 

= zfi, + , (122) 

the hrst term denoting the rotation about the precessing axis z , and the second term being 
the precession rate of z relative to the inertial frame (X, Y, Z) . This precession rate is 
given by 


^{^nert^al) ^ COS hp , jp siu hp , hp ) ^ , (123) 

because this expression satishes x z = z . 

In a frame co-precessing with the equator of date, the precession rate will be represented 
by vector 


= Ri_ 




(inertial) 


(124) 


where the matrix of rotation from the equator 
frame to the precessing one) is given by 

cos hp 
— cos Ip sin hp 
sin Ip sin hp 


of epoch to that of date (i.e., from the inertial 



sin hp 

0 

cos Ip 

sin hp 

sin Ip 

— sin/p 

sin hp 

cos/p 



From here we get the components of the precession rate, as seen in the co-precessing coordinate 
frame {x , y , z) : 

/X = (/ii , /Xs , h 3 ) = (jp 1 hpSin/p , hp cos Ip) . (125) 


28 

Tft€LLhjCJM.€E.£tc.€LL "Ph.y.A.Lc.A 


References 


[1] Berkovich, L. M. 2002. “Research of Non-Stationary Problems of Celestial Mechanics Trans¬ 
formation to an Autonomous Form.” In: Non-Stationary Dynamical Problems in Astronomy. 
Collection of articles, ed. by T. B. Omarov. Nova Science Publishers, New York, 2002. 

[2] Brumberg, V. A., L. S. Evdokimova, & N. G. Kochina. 1971. “Analytical Methods for the 
Orbits of Artihcial Satellites of the Moon.” Celestial Mechanics, Vol. 3, pp. 197 - 221. 

[3] Colombo, G. 1966. “Cassini’s second and third laws.” The Astronomical Journal, Vol. 71 
pp. 891 - 896. 

[4] Efroimsky, Michael. 2002a. “Equations for the orbital elements. Hidden symmetry.” 
Preprint 1844 of the Institute of Mathematics and its Applications, University of Minnesota 
http: //www. ima.umn.edu/preprints/feb02/feb02.html 

[5] Efroimsky, Michael. 2002b. “The Implicit Gauge Symmetry Emerging in the N-body Prob¬ 
lem of Celestial Mechanics.” astro-ph/0212245 

[6] Efroimsky, Michael, & Peter Goldreich. 2003. ’’Gauge Symmetry of the N-body Problem in 
the Hamilton-Jacobi Approach.” Journal of Mathematical Physics, Vol. 44 , pp. 5958 - 5977 
astro-ph/0305344 

[7] Efroimsky, Michael, & Peter Goldreich. 2004. ’’Gauge Freedom in the N-body Problem of 
Gelestial Mechanics.” Astronomy & Astrophysics, Vol. 415 , pp. 1187 - 1199 
astro-ph/0307130 

[8] Efroimsky, M. 2004. “Long-term evolution of orbits about a processing oblate planet. The 
case of uniform precession.” astro-ph/0408168 

This preprint is a very extended version of the published paper Efroimsky (2005). It contains 
all technical calculations omitted in the said publication. 

[9] Efroimsky, M. 2005a. “Long-term evolution of orbits about a processing oblate planet. 1. 
The case of uniform precession.” Celestial Mechanics and Dynamical Astronomy, Vol. 91 , 
pp. 75 - 108 

[10] Efroimsky, M. 2005b. “Long-term evolution of orbits about a processing oblate planet. 2. 
The case of variable precession.” In preparation. 

[11] Euler, L. 1748. Recherches sur la question des inegalites du mouvement de Saturne et de 
Jupiter, sujet propose pour le prix de I’annee. Berlin. Modern edition: L. Euler Opera 
mechanica et astronomica. Birkhauser-Verlag, Switzerland, 1999. 

[12] Euler, L. 1753. Theoria motus Lunae exhihens omnes ejus inaequalitates etc. Impensis 
Academiae Imperialis Scientarum Petropolitanae. St.Petersburg, Russia 1753. Modern 
edition: L. Euler Opera mechanica et astronomica. Birkhauser-Verlag, Switzerland 1999 

[13] Fukushima, T., & Ishizaki, H. 1994. “Elements of Spin Motion.” Celestial Mechanics and 
Dynamical Astronomy, Vol. 59 , pp. 149 - 159. 

[14] Goldreich, P. 1965. “Inclination of satellite orbits about an oblate processing planet.” The 
Astronomical Journal, Vol. 70 , p. 5 - 9 


Tft€LLhjc.§pL€E.£LC.€M.L 


29 



[15] Gurfil, P., and Klein, I. 2006. “Mitigating the Integration Error in Nnmerical Simnlations 
of Newtonian Systems.” Snbmitted to The International Journal for Numerical Methods in 
Engineering. 

[16] Gnrfil, P. 2004. “Analysis of J 2 ~pertnrbed Motion using Mean Non-Osculating Orbital 
Elements.” Celestial Mechanics & Dynamical Astronomy, Vol. 90 , pp. 289 - 306 

[17] Kinoshita, T. 1993. “Motion of the Orbital Plane of a Satellite due to a Secular Ghange 
of the Obliquity of its Mother Planet.” Celestial Mechanics and Dynamical Astronomy, Vol. 
57 , pp. 359 - 368 

[18] Lagrange, J.-L. 1778. Sur le Probleme de la determination des orbites des cometes d’apres 
trois observations, 1-er et 2-ieme memoires. Nouveaux Memoires de I’Academie de Berlin, 
1778. Later edition: Qluvres de Lagrange. Vol. IV, Gauthier-Villars, Paris 1869. 

[19] Lagrange, J.-L. 1783. Sur le Probleme de la determination des orbites des cometes d’apres 
trois observations, 3-ieme memoire. Ibidem, 1783. Later edition: Qluvres de Lagrange. 
Vol. IV, Gauthier-Villars, Paris 1869. 

[20] Lagrange, J.-L.1808. “Sur la theorie des variations des elements des planetes et en par- 
ticulier des variations des grands axes de leurs orbites”. Lu, le 22 aout 1808 a I’Institut de 
France. Later edition: Qluvres de Lagrange.Vol.Yl, pp.713-768, Gauthier-Villars, Paris 1877 

[21] Lagrange, J.-L. 1809. “Sur la theorie generale de la variation des constantes arbitraires 
dans tons les problemes de la mecanique”. Lu, le 13 mars 1809 a I’Institut de France. Later 
edition: (Euvres de Lagrange.Vol.Yl, pp.771-805, Gauthier-Villars, Paris 1877. 

[22] Lagrange, J.-L.1810. “Second memoire sur la theorie generale de la variation des constantes 
arbitraires dans tons les problemes de la mecanique”. Lu, le 19 fevrier 1810 a I’Institut de 
France. Later edition: Qluvres de Lagrange. Vol. VI, pp.809-816, Gauthier-Villars, Paris 1877 

[23] Lainey, V.; Gurhl, P.; & Efroimsky, M. 2005. “Long-term evolution of orbits about a 
processing oblate planet. 3. A semianalytical and a purely numerical approaches.” 

In preparation. 

[24] Mysen, E. 2004. “Rotational dynamics of subsolar sublimating triaxial comets,” Planetary 
and Space Science, Vol. 52 , pp. 897 - 907. 

[25] Newman, W., & M. Efroimsky. 2003. “The Method of Variation of Gonstants and Multiple 
Time Scales in Orbital Mechanics.” Chaos, Vol. 13, pp. 476 - 485. 

[26] Plummer, H. G. 1918. An Introductory Treatise on Dynamical Astronomy. Gambridge 
University Press, UK. 

[27] Slabinski, V. 2003. ’’Satellite Orbit Plane Perturbations Using an Efroimsky Gauge Ve¬ 
locity.” Talk at the 34th Meeting of the AAS Division on Dynamical Astronomy, Gornell 
University, May 2003. 

[28] Zanardi, M. G.; Vilhena de Moraes, R. 1999. “Analytical and semi-analytical analysis of 
an artihcial satellite’s rotational motion.” Celestial Mechanics and Dynamical Astronomy, 
Vol. 75 , pp. 227 - 250. 


T^uJmjc. Tft€LLhjc.§pL€E.£LC.€M.L 


30 



“Pu^UC. Tft€LLhjC.§M.€E.£LC.€M.L “Ptt.y.A.LC.A. 



